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Abstract 

A new research area of multiple methods integration is proposed for joining diverse methods 
of structural mechanics analysis which interact with one another. Three categories of multiple 
methods are defined: those in which a physical interface is well-defined; those in which a physical 
interface is not well-defined, but selected; and those in which the interface is a mathematical 
transformation. The integration procedures required in each of these categories is of necessity 
generally different. It is discussed how multiple methods integration is needed in order to effec- 
tively use analysis and design methods in the most appropriate structural region and to enable 
hierarchial analysis and design so that the effect of component design changes on full vehicle 
performance can be readily addressed. Two fundamental integration procedures are presented 
which can be expanded to integrate various methods ( e.g. t finite elements, Rayleigh Ritz, 
Galerkin and integral methods) with one another. Since the finite element method will likely be 
the major candidate method to be integrated, its enhanced robustness under element distortion 
is also examined and a new robust shell element is demonstrated. 

Introduction ^ 

Over the past three decades, finite elements have become the workhorse of structural analysis. 
This is due in part to the ability of finite elements to allow modeling of arbitrary and complex 
geometric shapes. With the increased utilization of composite aerospace structures, there is an 
even greater desire for modeling flexibility. Manufacturing of composite aerospace components, 
leading to new geometric configurations, and the need for detailed modeling as required for 
composite failure prediction are both producing a desire for added modeling flexibility. 

Though finite elements dominate the structural analysis discipline, this does not imply that 
other modeling methods are not valuable. To the contrary, special purpose modeling techniques 
and their associated software abound throughout the industry because of their proven utility. 
However, little has been done to combine these modeling methods with one another and with 
finite element modeling. Consequently, the special purpose analyses are rarely used in a hierarchial 
sense where analyses at different levels of detail or on different structural component levels inherit 
data from the other levels and interact with one another. One purpose of this paper is to present 
a relatively new endeavor referred to as multiple methods integration which will enable hierarchial 
analysis and to present some illustrative examples with fundamental equations for the integration 
process. Giles and Norwood[l] initiated this area of research thus, the intent here is not to claim 
a new concept, but to build on their work and to continue to motivate research which will lead 
to practical application. 

Since the finite element method will likely be one of the multiple methods which is integrated 
with other methods, it is desirable to examine the finite element capability in the light of its 
readiness for integration since an integrated method will be no better than its constituents. Thus, 
a second purpose of this paper is to examine the robustness of the finite elements themselves 
which provides the modeling flexibility of the finite element method. In spite of the plethora 
of finite elements types presently available, (see references [2] and [3] for a broad review), the 
need still exists for robust elements which provide adequate performance even when the finite 
elements have distorted shapes. This paper reviews the distortion resistant capability of several 
types of shell finite elements as well as a newly developed shell element which is very distortion 
insensitive. 
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Multiple Methods Integration 

Over the years, many special purpose structural analysis methods have been developed for 
various classes of structural components under a variety of loading conditions. Also, the finite 
element method has emerged as the leading method for general purpose analysis. Just as the 
various components comprising an aerospace vehicle must be integrated together, so must the 
methods developed for their analysis. Herein this process is referred to as multiple methods 
integration. For several reasons it is advantageous to be able to readily and efficiently merge 
or integrate the various special purpose analysis methods with each other and with the general 
purpose finite element method. Among the reasons for multiple methods integration are the 
following: utilization of the preferred analysis method in the most appropriate structural region; 
utilization of high performance computers which operate in parallel on different regions and 
on different methods; retention of modularity in software architecture; refinement of structural 
regions independently of one another; and hierarchial analysis. It is expected that multiple 
methods integration should lead to more accurate and efficient analysis and design, since the 
total analysis becomes more tailored to the article being analyzed. Even more important than 
reduced computer time, however, is the promise that multiple methods integration can reduce 
the time spent by the analyst or designer. 

If done properly, the integration procedure will allow the designer to modify the design and/or 
refine the model of one structural article independently of models of other articles. This is 
not now generally available. The effect of a redesigned component may propagate to other 
components as well as to the full vehicle. Predicting this propagation and the trade-offs it 
causes is critical to efficient design. However, not only does the redesign propagate to the rest 
of the vehicle, but in general, so does the remodeling which is done to account for the redesign. 
It is difficult to separate these two effects. Multiple methods integration should be such so as 
to allow the remodeling effects to be filtered out of the redesign. 

On the other hand, certain affects may be very localized, and in accordance with St. Venant’s 
principle do not propagate far. These can be handled by specialized methods and then interfaced 
to methods treating the global response. 

Once the integration procedures are established, the COMET (computational MEchanics 
Testbed) at LaRC (reference [4]) can be used as the vehicle for implementing the integration 
and validating the integration procedure. The special purpose computer software associated 
with each method can then be readily integrated with one another and with the finite element 
software capability residing in COMET. 

In this section three categories of multiple methods are outlined and some concepts for integrating 
the methods are presented. These categories are identified by the type of interface existing 
between the multiple methods. They are as follows: (1) the case where two or more methods 
are employed in different structural regions and the physical interface between the regions is 
known or well-defined; (2) simlar to the first case in that the interface is physical, but it cannot 
be defined, is not known or is selected by the analyst to produce an efficient model; and (3) 
the interface is not physical, but mathematical. Examples of the three categories are given in 
Figure 1. 

Category 1. Physical Interface Well-Defined - In this category there are several well-known pro- 
cedures for integration, such as, component substructuring, modal synthesis and interdisciplinary 
analysis such as fluid-structure analysis. Except in the case of interdisciplinary analysis, the gen- 
eral practice is to use the same methods in each structural region. Little has been done to 
interface different types of methods such as finite elements with Galerkin methods. Moreover, 
modal synthesis methods still have serious shortcomings and little has been done in the non-linear 
area for this entire category of methods. 

Category 2. Physical Interface Not Well-Defined - In this category as in the first, the interface 
is physical; however, it is not the intersection between two structural components. Rather, the 
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interface Is chosen by the analyst In order to produce a more efficient model of a single structural 
component. For example, this category includes such techniques as, global-local modeling where 
there are near and far field models, multiple transient analysis modeling where different time steps 
or different temporal Integrators may be used in different regions, and 2D-to-3D modeling where 
plate finite elements may be tied to brick finite elements. Often, the accuracy of the model is 
highly dependent on the placement of the interface. 

Category 3. Mathematical Interface - In this category two or more methods are used to model 
the same structural article and hence there is no physical interface. Rather, the interface is 
mathematical and often represents the transformation(s) between two or more modeling phases 
of the analysis. Examples in this category are methods such as predictor/corrector, hybrid, re- 
analysis multi-grid methods and reduced basis methods. An example of a category 3 multiple 
method is depicted in Figure 2. References [5-8] develop this procedure wherein finite element 
generated results are enhanced by a process called re-continualization. This process derives a 
functional fit to the discrete finite element results. Once the functional fit is derived, detail 
stresses may be calculated by differentiating the function. The critical aspect of the process in 
deriving the function is that no differencing of the discrete finite element results is required. 




(b) Category 2 
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Figure 1. Three Categories of Multiple Methods. 
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Multiple methods integration enables hierarchial analysis by tying together the results of various 
levels of structural detail. This can be accomplished in category (1) multiple methods where 
results from one component level become input to the next level and in category (3) multiple 
methods wherein levels of interacting methods treat levels of progressively complex structural 
detail. This is especially needed in composites since detail stresses are required for failure 
prediction. 
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Figure 2. New Method Developed for Enchancing Finite Element Derived Results. 

The remainder of this section concentrates on two integration methods for categories (1) and 
(2). Whereas, the two concepts laid down here are applicable to both categories, one tends to be 
more applicable to category (1) and the other to category (2). In category (1) multiple methods 
it is often critical to model the interface with sufficient refinement to predict accurate stresses, 
since the interface between two components is a stress critical region where failures usually 
initiate. The model of the interface for category (1) tends to refine as the design matures. On 
the other hand, category (2) multiple methods do not necessarily seek accurate stresses at the 
interface since they are used where the stress field does not contain steep stress gradients. The 
category (2) integration should be such that the results are insensitive to the precise placement 
of the interface. 

The concept of each method is described by considering two structural regions having a common 
interface (see Figure 3). In these concepts, the displacements and/or normal gradients along 
the interface are assumed to be described by some functional forms. The two concepts are 
distinguished by the way the models for the two regions are tied to the functional forms. In the 
first concept, the models are rigidly tied to the functional form whereas in the second they are 
loosely tied to the functional form by a least squares procedure. As stated in the introduction, 
the intent here is not to claim some new concepts, but rather to stimulate research into the large 
variety of integration procedures possible so as to develop methods for practical application. 
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Applied Displacements 



Figure 3. Two Connected Regions. 


Method 1: Rigid Connection via Functional Form 

Consider an interface between two structural regions. A displacement or normal gradient com- 
ponent along this interface is assumed to have a functional form given by 


q = q(s) (1) 

where s is the interfacial path coordinate which is defined in terms of the analysis coordinate 
system. The functional form for q(s) uses unknown coefficients and when the function is evaluated 
at a discrete coordinate, one has, 

q(si) = Ua ( 2 ) 

where, for a polynomial of order M - 1, 


U 


(1, 


> °i * * ' 




M- 




(3) 


and a is a column vector of coefficients for the polynomial q(s). 

It is assumed that the values of 9 associated with the model of region (1) and those of region 
(2) can be expressed as two sets of discrete points lying along the common interface between 
the regions. These discrete values are then slaved to the assumed functional form as, 


q(*>)-q i 1} = 0 and «(«>)- 9^ = 0; s = 1,2 IV 


«(*)-, 




3 = 1,2 !V< a) 


(4) 


where, 9,^ and q are discrete values of q associated with region (1) or (2) respectively, at 
values of and Sj. (The values of Si and Sj need not be coincident.) Then, 


<?<*) = T(*>o; *=1,2 (5) 

where, is a column vector associated with the interface of region (*), all of whose entries 
are of one type of degree-of-freedom (either a displacement component or its gradient). The 
matrix Q[ k ^ is 1 x 1 and is to be made continuous across the interface. The matrix is 

N < * ) x M and its rows are U, Equations (5) represents constraints and the functional form 
introduces M new unknown variables, thus there are - M interface constraints using 

this interfacing technique. 
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The order of the polynomials used to represent different types of degrees-of-freedom need not 
be the same, but for simplicity, they are assumed to be of the same order in this paper. Hence, 
for all degrees-of-freedom which are to be made continuous across the interface, 


qW- T (k) a . k = l,2 ( 6 ) 

where Q ^ is an extended column vector of all degrees-of-freedom which are to be continuous 
across the interface and T is a block diagonal matrix whose blocks are T* . Thus, Q ^ is 
jy(*) x 1, T< fc) is NW x M and a is M x 1 where N W is nlf k) with n the number of different types 
of degrees-of-freedom which are to be made continuous across the interface. 

Variational Equations - The models associated with each of the regions yields the virtual work 
of each as, 

= 6D (k) K (k) D {k) ( 7 ) 

where is the column vector of all degrees-of-freedom on the boundary of region (k) and 
is the associated stiffness matrix for region (fc) which is typically derived from the full stiffness 
matrix by reducing out the interior degrees-of-freedom. It could be generated by a finite element 
model or some other model. The vector £)(*) is now partitioned into three parts; those N free 
degrees-of-freedom lying on the interface of region (fe); those L free degrees-of-freedom lying 
elsewhere on the boundary of region ( k ), and those P<*> degrees-of-freedom on the boundary 
of region (ifc) which correspond to fixed or applied boundary displacements. After partitioning, 
equation (7) becomes, 
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( 8 ) 


For the case of no applied external forces, 

6UW + 6UW = 0 ( 9 ) 

Since equation (9) must be valid for all variations, the governing equation for the case of a rigid 
interface connection via a functional form is, 
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where, 

Sn = t (i > t 4/ )t(1 ) + t< j)T 4/ )t(2) 

S/1 = 

S n = 

S n = 4 f 

Sjj — Kpp 

Solution of equation (10) for a, and Q ^ gives, 

a = S- 1 [(S/xS- 1 ^] - k< ia)Qa ] + (Sn&K™ - K%)q { ?} 

Qf ] = -5u l [k'fW? + Sfra] 

Q^ = -S^[K^ A Q^ + Sj 2 a] 


(11) 


where, 


S = S n - Sn $ii Sn - SnS^Sf, 
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Thus, the procedure for a rigid interface connection via a functional form is 

1. Model regions 1 and 2 independently 

2. Create stiffness matrix with all degrees-of-freedom lying interior to the boundaries reduced 
out 

3. Calculate a, and Q {2) from equations (11) 

4. Recover interface values from equations (6) 


Method 2: Least Squares Interface Connection via a Functional Form 

In this method, equation (4) is replaced by the squared equations, 

n< 1) 2 . 

£[«(*> -«!'’] =«and £ [,(,) -,}’»] =0 (12) 

» = 1 j = 1 

Using the process of least squares and equation (6) produces 2Af equations, namely, 

T( 1 ) T (T< 1 >a - Q ( j >) = 0 and (T (J >a - Q ( p) = 0 
which yields M constraint equations between and \ namely, 

^7 , ( 1 ) t 7’( 1 )^ j , ( 1 ) t q( 1) _ ^y(a) T T(a)J" l | r (a) T g(> ) = 0 


or, 

CD = 0 (13) 

where, C is a constraint matrix of order M x ( N W 4- 4 L 4 i( 2 )) expressed as, 


c= [r^o.-r^o]; 


d = 
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Returning to equation (7) and applying the constraints via Lagrangian multiplier technique yields, 
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(14) 


subject to the constraints, 


CD = 0 

With A a Lagrangian multiplier vector, the term C T X in equation (14) represents the interfacial 
forces between regions (1) and (2). From equation (13), these forces are orthogonal to the 
displacement vector D which does not violate the constraints. 

Method (2) produces one interface constraint for each polynomial coefficient. That is, the 
number of constraints is Af. This is in contrast with method (1) which produces 4 - M 

constraints. 
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Equation (14) may be solved by transforming to a reduced set of response variables as, 


D-BD (15) 

with B of the order N x (N - R) where N is N M + N and R is the rank of the 
matrix E = C T C. If all the constraints are linearly independent, than the rank R will be equal to 
the number of constraints, namely, M. Under certain circumstances, the constraints may not be 
all independent. For example, if some gridpoints along the interface are coincident and are also 
rigidly constrained to one another. This difficulty is handled by establishing B through a singular 
value decomposition (SVD) which produces B such that its columns are the eigenvectors of E 
associated with zero eigenvalues. Substituting equation (15) into equation (14), pre-multiplying 
by B t and solving for D yields 

D = (B T kB)~ 1 B T Q (16) 

where K and Q are the matrix and column vectors of equation (14) and the Lagrangian mul- 
tiplier term has vanished since B T C T vanishes. Equation (15) is then used to recover the full 
displacement vector D. 

It is not difficult to extend the integration methods (1) and (2) to the case where the analysis 
methods of the regions are of the global Ritz or Galerkin type. It is quite likely though, that one 
region will be modeled using finite elements especially if this region is characterized by complex 
geometry, for the finite element method is typically thought of as having the positive attribute of 
treating arbitrary geometry. However, arbitrary geometry often leads to distorted finite elements. 
Unfortunately most finite elements are not robust to element distortion, (/.e., they are distortion 
sensitive). Thus attention is next focused on the robustness issue since it will be a critical aspect 
of making multiple methods integration a practical reality. 

Robustness of Finite Elements 

The finite element modeling of general aerospace shell structures usually requires the use of 
distorted meshes. Additional complexities arise if the original finite element discretization is 
further refined through an adaptive refinement procedure. Most of the 4-node shell elements 
developed in the past do not produce reliable results for distorted meshes. Whereas, the element 
developers are usually aware of the limitations and pitfalls of the elements, the users may not 
be aware of all these limitations and may make invalid use of them. Therefore, it is desirable to 
formulate simple 3-node and 4-node shell elements that are free from the usual limitations and 
pitfalls, such as locking, sensitivity to mesh distortion, non-invariance, and spurious modes, and 
more importantly produce accurate and reliable results. 

One way of attacking the shortcomings of membrane elements is to include the nodal normal 
rotational or “drilling" degrees of freedom in the element formulation. In early attempts, these 
rotational degrees of freedom were used in cubic displacement functions. However, Irons and 
Ahmad demonstrated that this approach had serious deficiencies^]. The elements formed in this 
manner force the shearing strain to be zero at the nodes, and because these elements do not 
pass the patch test, they could produce erroneous results in some structural analysis problems. 
Recently researchers have used these rotational degrees of freedom in quadratic displacement 
functions with more success[10-15]. In previous papers, this latter method has been employed 
in the following way. First, the element is internally assumed to be an 8-node isoparametric 
element with 4 corner nodes and 4 midside nodes each having two displacement degrees of 
freedom, and the stiffness matrix associated with this “internal'' element is calculated. Then, 
this stiffness matrix is condensed to that corresponding to a 4-node element with 12 degrees 
of freedom by associating the displacement degrees of freedom at the midside nodes with the 
displacement and rotational degrees of freedom at the corner nodes. MacNeal[12] has used this 
approach to develop a 4-node displacement-based membrane element with selective reduced- 
order integration. Yunus et al., in reference [13], have also used this method to develop an 
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assumed-stress hybrid/mixed membrane element. Aminpour[14] has also used this method to 
develop an assumed-stress hybrid/mixed shelf element. 

In this paper, a 4-node assumed-stress hybrid quadrilateral shell element with “drilling” degrees 
of freedom is presented. The formulation is based directly on a 4-node element from the 
outset in contrast to elements whose formulations began with an “internal” 8-node element. 
Formulating the element in this manner bypasses the formation of the stiffness matrix for an 
8-node isoparametric element and the subsequent transformation of this stiffness matrix to that 
corresponding to the stiffness matrix of a 4-node element. 

Hybrid Variational Principle 

The classical assumed-stress hybrid formulation of Pian[16] is based on the principle of minimum 
complementary energy. The displacements are described on the element boundary and an equili- 
brating stress field is described over the the domain of the element. It was later recognized that 
the same method may be derived from the Hellinger-Reissner principle ( e.g. f see reference [17]). 
However, in the Hellinger-Reissner principle, the stress field does not have to satisfy the equilib- 
rium equations a priori, and the displacement field has to be described over the domain of the 
element and not just on the boundaries. The membrane element in reference [13] and the shell 
element in reference [14] were both developed using the Hellinger-Reissner principle. However, 
an assumed-stress hybrid 4-node shell element similar to that of reference [14] may also be easily 
formulated using the minimum complementary energy principle with the advantage being that 
only the displacements on the boundary of the element enter into the formulation. As such, the 
formulation is then based directly on a 4-node element rather than internally formulated as an 
8-node element and then condensed to a 4-node element. 

The invariant properties of the element are preserved by the proper choice of a local element 
Cartesian-coordinate system. The local element Cartesian-coordinate system is shown in Figure 4 
and is obtained by bisecting the diagonals of the element. The axes of this coordinate system are 
approximately parallel to the edges of the element for non-rectangular geometries (e.g., tapered 
and skewed elements) which would make the element less sensitive to mesh distortion. 

The assumed-stress hybrid formulation is based on assuming an equilibrating stress field <x in the 
interior of the element as a = Pf. 3, and assuming the displacement field u only on the boundary of 
the element as u = Nq. The matrices P and N consist of the interpolating functions for stresses 
and displacements, respectively, and coefficients (3 and q are the unknown stress parameters and 
nodal displacements and rotations, respectively. 



Figure 4. Element coordinate system definition. 

Displacement Field Description 

All three displacement components on the element boundary are assumed to vary quadratically 
and all three rotational components to vary linearly. The expressions for these boundary dis- 
placements and rotations were derived in detail in reference [14] and only the final results are 
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given herein. The in-plane boundary displacements on edge 1 of Figure 4 are given by 

u = -(1 - £ )i*i 4- -(1 + -+ — jp(l - - M 

228 ( 17 ) 

v = ^(1 - £)vi + ^(1 + {) v 2 - — jp(! - - **l) 

and the out-of-plane boundary displacement and rotations on edge 1 of Figure 4 are given by 
w = i(l - £)wi + ^(1 +0 U '2 jp(l _ £ 3 )(^*2 _ ^*0 + _ ^yi) 

** = 5 ( 1 - 0 **i + ^(i + O 0 *j ( 18 ) 

fa = 2 ( 1 - O^yi + 2 ^ + ^y J 

where, A*i and Ayi are the Ax and Ay of edge 1 with respect to the reference local element 
x-y coordinate system (e.g., A*i = x 2 — * 1 ) and £ is a non-dimensional coordinate on edge 1 
such that £ = -1 at node 1 and £ = +1 at node 2. The description for the displacements and 
rotations on the other edges of the element are readily obtained. 

The displacement and rotation descriptions in equations (17) and (18) allow for in-plane shearing 
strain and transverse shearing strains, respectively. This feature is in contrast to cubic interpo- 
lations of in-plane or out-of-plane displacements which force the in-plane shearing strain or the 
transverse shearing strains to be zero at the element nodes. The elements using cubic interpo- 
lation do not pass the patch test and perform poorly for some structural analysis problems[9]. 
The elements constructed using quadratic interpolation, on the other hand, pass the patch test 
which is a necessary condition for convergence to the correct solution. 

Stress Field Description 

The stress field should be selected in such a manner that no spurious zero-energy mode is 
produced. In order to avoid spurious zero-energy modes, each independent stress term must 
suppress one independent deformation mode. Therefore, the minimum number of stress terms 
required is equal to the number of degrees of freedom of the element less the number of 
rigid body modes. The equilibrating stress field is expressed in the local element Cartesian- 
coordinate system shown in Figure 4 and is similar to that proposed by reference [14]. However, 
in reference [14] the Hellinger-Reissner principle was used, and the stresses were expressed in the 
natural-coordinate system. 

The following equilibrating stress (resultant) field is considered for the membrane part 

N X =P i+fay + fax + fay 2 

N y = fa +fax + fay + fax 2 ( 19 ) 

N xy = fa- fay - fax 

The first five terms of the stress field in equation (19) represent the stress field that was used in 
the original 4-node (see reference [16]) assumed-stress hybrid membrane element with 8 degrees 
of freedom which did not include any normal rotational degrees of freedom. The remaining four 
terms are present to suppress the four rotational degrees of freedom present in this formulation. 

The following equilibrating stress (resultant) field is selected here for the bending part 

M x = fa + fay + fax + faxy 

My = fa + fax + fay + faxy 

M zy — fa + faox + fax y+ -Aa* 2 + -/?i3y 3 

Qx = {fa + /3n) + (fa + Pi3)y 
Qy — (fa "I" Ao) + (^e + fai) x 
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The stress fields in equations (19) and (20) produce no spurious zero-energy modes[14] and [15]. 
It is observed that both the membrane and bending stress (resultant) fields remain invariant upon 
node renumbering. 


NUMERICAL RESULTS 

The performance of the 4-node quadrilateral shell element developed in this paper (herein referred 
to as AQD4) is evaluated in this section. The element has been implemented in the NASA 
Langley COMET software system[4]. Selected test problems in this article are the straight 
cantilever beam and the Scordelis-Lo roof. More extensive test problems for this element are 
reported in reference [15]. The results for the present element are compared with the results 
using the QUAD4 element of the MSC/NASTRAN from reference [20], the Q4S element from 
reference [12], the ES1/EX47[18], and ES5/E410[19] elements of the NASA Langley CSM 
Testbed. A brief description of these elements is given in reference [15]. The dimensions and 
properties for the test problems are chosen in consistent units. 

Straight Cantilever Beam 

The straight cantilever beam problem suggested in reference [20] is solved for the three dis- 
cretizations shown in Figure 5. Normalized results for the present element along with the results 
for other elements are shown in Table 1. These results indicate that all elements perform well 
for the rectangular mesh. However, for the trapezoidal and parallelogram meshes which contain 
considerable amount of distortion, only the Q4S and the present element (AQD4) perform well. 
The present element produces an error of less than 3.5% for all meshes and loads which indicates 
the insensitivity of the present element to mesh distortion. 
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Figure 5. Straight cantilever beam problem. Length=6., height=0.2, depth=0.1, E=10 7 , u-0.3, 
mesh=6xl. Loading: unit forces at the free end. 
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Table 1. Normalized tip displacements in direction of 
loads for straight cantilever beam. 


Tip Loading 
Direction 

QUAD4 

MSC/ 

NASTRAN 


ES5/ 

E410 

(STAGS)t 

AQD4 

(present) 

(a) rectangular shape elements 

Extension 

.995 

.995 

.994 

.998 

In-plane Shear 

.904* 

.904 

.915 

.993 

Out-of-Plane Shear 

.986 

.980 

.986 

.981 

Twist 

.941** 

.856 

.680 

1.011 

(b) trapezoidal shape elements 

Extension 

.996 

.761 

.991 

.998 

In-plane Shear 

.071* 

.305 

.813 

.986 

Out-of-Plane Shear 

.968 

.763 

# 

.965 

Twist 

.951** 

.843 

# 

1.009 

(c) parallelogram shape elements 

Extension 

.996 

.966 

.989 

.998 

In-plane Shear 

.080* 

.324 

.794 

.972 

Out-of-Plane Shear 

.977 

.939 

.991 

.980 

Twist 

945 ** 

.798 

.677 

1.010 


f These elements are not invariant and do not pass the patch test. 

* The results from MacNeal’s new Q4S element (not shown in table) for in-plane shear load are 
reported in reference [12] to be .993, .988, and .986 for the meshes (a), (b), and (c) in Figure 5 
respectively. 

** These results for twist were normalized with .03028 in reference [20]. Herein, all the other 
results for twist are normalized using .03046 according to Timoshenko and Goodier's Theory of 
Elasticity. 

# The element produces a singular stiffness matrix for this mesh. 


Scordelis-Lo Roof 

The Scordelis-Lo roof is shown in Figure 6 . This structure is a singly-curved shell problem, 
loaded under gravity, in which both the membrane and bending contributions to the deformation 
are significant. Because of symmetry, only one quadrant of the problem is modeled. The mesh 
on one quadrant is chosen to be NxN for N=2, 4,6,8,10 (N=number of elements along each 
edge) to show the convergence of the solutions for the present element. Normalized results are 
shown in Table 2. For this problem, the mesh is made of uniform rectangular-shaped elements 
and all the elements in the table perform well. It is observed that the convergence rate to the 
reference solution for the present element is roughly the same as the other elements and thus, 
the addition of the rotational degrees of freedom does not affect the convergence rate of the 
present element for this problem. 
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Figure 6. Scordelis-Lo roof problem. Length=50., radius=25., thickness=0.25, E=4.32xl0 8 , 
0.3, mesh=NxN. Loading: 90. per unit area in vertical direction, /.e. f gravity load; 
u x =ti z =0 on curved edges. Reference solution: vertical displacement at midpoint of 
free-edge=0.3024 from reference [20]. 


Table 2. Normalized displacements at the midpoint of 
the free-edge for Scordelis-Lo roof. 


Mesh 

QUAD4 

MSC/ 

NASTRAN 

ESI 

EX47/ 

(ANS) 

ES5 

E410/ 

(STAGS) 

AQD4 

(present) 

2x2 

1.376 

1.387 

1.384 

1.218 

4x4 

1.050 

1.039 

1.049 

1.021 

6x6 

1.018 

1.011 

1.015 

1.006 

8x8 

1.008 

1.005 

1.005 

1.003 

10x10 

1.004 

1.003 

1.001 

1.001 


Concluding Remarks 

A new research area of multiple methods integration is proposed for joining diverse methods 
of structural mechanics analysis which interact with one another. Three categories of multiple 
methods are defined: those in which a physical interface is well-defined; those in which a physical 
interface is not well-defined, but selected; and those in which the interface is a mathematical 
transformation. The integration procedures required in each of these categories is of necessity 
generally different. 

Two fundamental integration procedures have been formulated. One which establishes a tight 
interface connection and finds application in the first category, and one which establishes a more 
forgiving connection and finds application in the second category. 

Moreover, since the most likely candidate for integration with other methods is the general 
purpose finite element method, the robustness issue of finite elements to element distortion has 
been examined and a new robust 4-node quadrilateral shell element demonstrated. 

The membrane part of the element formulation includes the drilling degrees of freedom and the 
bending part is of class C° and takes into account the effects of transverse shear deformations. 
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The element formulation is simple and straightforward, and is derived directly for a 4-node 
element rather than being reduced from “internal" 8-node isoparametric elements. Thus the 
formation of the element kernel does not require the time consuming transformation from an 
8-node to a 4-node element. 
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